Jump to Empirical Spatial Indicator Analysis

Jump to Location
Jump to Range
Jump to Occupancy
Jump to Aggregation

1 Introduction

This vignette provides an overview of spatial indicators used to estimate the location, range, occupancy, and aggregation of populations in the context of fisheries management.

1.1 Data Preparation

Spatial indicators are configured to work with DATRAS survey data. Load in some survey data and process it first before parsing through spatial indicator functions.

Load DATRAS exchange data

Here we load NS-IBTS Quarter 1 survey data for years 2018 to 2022.

yrs <- c(2000:2022)
qrs <- c(1)
srv <- "NS-IBTS"

hh  <- icesDatras::getDATRAS(record = "HH", srv, years = yrs, quarters = qrs)
hl <-  icesDatras::getDATRAS(record = "HL", srv, years = yrs, quarters = qrs)
ca  <- icesDatras::getDATRAS(record = "CA", srv, years = yrs, quarters = qrs) # for completeness, but not used

Remove duplicate rows in exchange data.

hh <- unique(hh)
hl <- unique(hl)
ca <- unique(ca)

ICES Statistical Rectangles

In the exchange data (in hh and hl) we know the ICES rectangle that each haul was conducted in, but we do not have the associated division. To be able to subset data to a specific stock region, we need to download ICES Statistical Rectangle shapefile (Quick Downloads > ICES StatRec mapped to ICES Area) and append information on ICES divisions to our exchange data. A download is available on GitHub and has been converted to a .rds file which is loaded in here.

load(url("https://github.com/peterjosephkidd/SpatIndAssess/raw/main/Data/ICES%20Rect/ices_rect.rds"))

Add ICES Divisions

area_div <- dplyr::distinct(ices_rect[c("ICESNAME", "Area_27", "Shape_Area")])
hh <- merge.data.frame(hh, area_div, by.x = "StatRec", by.y = "ICESNAME")
ca <- merge.data.frame(ca, area_div, by.x = "AreaCode", by.y = "ICESNAME")

Basic Data Processing

Remove invalid hauls:

hh <- filter(hh, !HaulVal %in% c("I", "P")) # p = partly valid, it is deprecated

-9 is a placeholder for NAs. Change to NA for TotalNo column:

hl$TotalNo[hl$TotalNo == -9] <- NA

Create haul.id:

hh$haul.id <- as.character(
  paste(hh$Year, hh$Quarter, hh$Country, hh$Ship, hh$Gear, hh$StNo, hh$HaulNo, 
        sep = ":"))
hl$haul.id <- as.character(
  paste(hl$Year, hl$Quarter, hl$Country, hl$Ship, hl$Gear, hl$StNo, hl$HaulNo, 
        sep = ":"))
ca$haul.id <- as.character(
  paste(ca$Year, ca$Quarter, ca$Country, ca$Ship, ca$Gear, ca$StNo, ca$HaulNo, 
        sep = ":"))

Merge haul information from hh with landings information in hl.

## Refine columns
m <- hh[c("haul.id", "Year", "Quarter", "Month", "Survey","Country", "Ship", 
          "Gear", "GearEx", "DoorType", "HaulDur", "HaulNo", "StNo", "SweepLngt", 
          "StatRec", "Area_27", "ShootLong", "ShootLat", "HaulVal", "Depth")]
## Merge HL with HH 
hlhh <- merge(hl, dplyr::distinct(m),
              c("haul.id", "Year", "Quarter", "HaulNo", "StNo", "Gear", "GearEx", 
                "DoorType","Ship", "SweepLngt", "Country", "Survey"))
rm(m)

1.2 Stock ple.27.420

Now that DATRAS survey data has been loaded and formatted, choose a stock to investigate and filter data accordingly. This example looks at European plaice (Plueronectes platessa) in Subarea 4 (North Sea) and Subdivision 20 (Skagerrak; ple.27.420).

Set parameters

# ple.27.420
species <- c("plaice", "Pleuronectes platessa")
species_aphia <- findAphia(species[2], latin = T)
stk_divs <- c("4.a", "4.b", "4.c", "3.a.20")

Filter data

Basic checks

2 Empirical Spatial Indicator Analysis

# Source Spatial Indicator functions from GitHub
i <- 1
source("https://raw.githubusercontent.com/peterjosephkidd/SpatIndAssess/main/Functions/spatinds_funs.R", print.eval = F)

2.1 Location

2.1.1 Centre of Gravity (CGx, CGy):

The mean longitudinal and latitudinal location of the population. Changes in CG indicate whether the mean location of the population is shifting eastward/westward (CGx) or northward/southward (CGy).

cg <- coginis(hlhh, yrs, qrs, species_aphia, stk_divs, 
              iso = F, inertia = F, # toggle off outputting two another spatial indicators for now
              density = T) # weight hauls by  density of catch at each sample site (F therefore uses binary presence-absence data)

head(cg)

  Year Quarter CoG (x) CoG (y)
1 2000       1    6.35    55.8
2 2001       1    4.74    54.7
3 2002       1    4.25    54.7
4 2003       1    4.16    54.6
5 2004       1    4.95    55.1
6 2005       1    5.87    55.8

Visualise changes in cg over time

cg_p <-cogplot(cg, 
               # set grid = ices_rect for ICES rectangles -- takes time to load
               grid = ices_rect, areas = stk_divs,
               xlim = c(-2, 6), ylim = c(52, 60)) 
suppressWarnings(ggplotly(cg_p, tooltip = "text"))

2.2 Range

2.2.1 Inertia (I):

Inertia describes the dispersion/variance of the population around its centre of gravity. High values of inertia indicate that the population is widely spread across space.

inert <- coginis(hlhh, yrs, qrs, species_aphia, stk_divs, 
              iso = F, inertia = T, # inertia toggled on
              density = T) #

head(inert)

  Year Quarter CoG (x) CoG (y) Inertia
1 2000       1    6.35    55.8    14.8
2 2001       1    4.74    54.7    12.3
3 2002       1    4.25    54.7    12.8
4 2003       1    4.16    54.6    13.5
5 2004       1    4.95    55.1    19.7
6 2005       1    5.87    55.8    19.6

Plot the trend over time

in_p <- ggplot(data = inert, aes(x = Year, y = Inertia)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Inertia (I) Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA))

in_p

2.2.2 Extent of Occurrence (EOO):

A convex hull polygon is drawn around occurrence points (i.e. sample sites with catch of the target species > 0). EOO is the area of the polygon. High EOO indicates that the population is spread over a large geographical area.

eoo <- chullarea(hlhh, yrs, qrs, species_aphia, stk_divs)[1:3] # this function needs tidying up, ignore other columns for now
head(eoo)

  Year Quarter convex_hull_area
1 2000       1             89.3
2 2001       1             85.9
3 2002       1             90.0
4 2003       1             87.0
5 2004       1             87.5
6 2005       1             88.3

Plot the trend over time

eoo_p <- ggplot(data = eoo, aes(x = Year, y = convex_hull_area)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  ylab("Extent of Occurence") +
  labs(title = "Extent of Occurence (EOO)",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA))

eoo_p

2.2.3 Ellipse Area (ELA):

Similar to EOO. But instead of a convex hull, ELA calculates the area of an ellipse that encompasses 95% of the occurrence sites.

print(stk_divs)

[1] "4.a"    "4.b"    "4.c"    "3.a.20"

ela <- ellarea(hlhh, yrs, qrs, species_aphia, stk_divs = stk_divs)
head(ela)

  Year Quarter Ellipse Area
1 2000       1         86.2
2 2001       1         75.0
3 2002       1         86.6
4 2003       1         92.0
5 2004       1        102.0
6 2005       1        109.4

ela_p <- ggplot(data = ela, aes(x = Year, y = `Ellipse Area`)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Ellipse Area (ELA) Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA))

ela_p

2.3 Occupancy

2.3.1 Proportion of Presence (POP):

Indicates the proportion of area occupied by the target species using binary presence-absence data.

2.3.1.1 Rectangle (POPR):

The proportion of the ICES rectangles surveyed with a catch of the target species > 0.

popr <- pa_rect(hlhh, yrs, qrs, species_aphia, stk_divs)
head(popr)

# A tibble: 6 × 5
   Year nrects Quarter nrects_p PosAreaR
  <dbl>  <int> <chr>      <dbl>    <dbl>
1  2000    163 1            127    0.779
2  2001    165 1            128    0.776
3  2002    162 1            131    0.809
4  2003    165 1            123    0.745
5  2004    165 1            131    0.794
6  2005    166 1            131    0.789

popr_p <- ggplot(data = popr, aes(x = Year, y = `PosAreaR`)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Proportion of Presence in ICES Rectangles (POPR) Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA)) +
  ylab("Proportion of Presence (Rectangle)")

popr_p

2.3.1.2 Haul (POPH):

The proportion of hauls with catch of the target species > 0.

poph <- pa_haul(hlhh, yrs, qrs, species_aphia, stk_divs)
head(poph)

# A tibble: 6 × 5
   Year no_haul.ids Quarter pr_hauls PosAreaH
  <dbl>       <int> <chr>      <dbl>    <dbl>
1  2000         367 1            256    0.698
2  2001         411 1            288    0.701
3  2002         401 1            297    0.741
4  2003         398 1            276    0.693
5  2004         355 1            246    0.693
6  2005         371 1            261    0.704

poph_p <- ggplot(data = poph, aes(x = Year, y = `PosAreaH`)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Proportion of Presence in Survey Hauls (POPR) Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA)) +
  ylab("Proportion of Presence (Haul)")

poph_p

2.4 Map Indicators

The mapdis() function visualises the above range and occupancy indicators on a single map for one timepoint (e.g. year and quarter). The function can equally be used to visualise a single indicator.

mapdis(hlhh, yrs = 2022, qrs, species_aphia, stk_divs, ices_rect, # data & parameters
       cog = T, inertia = T, EOO = T, ELA = T, # indicator toggles
       density = T, # weight samples 
       title = "Plaice (Pleuronected platessa)\nNS-IBTS", 
       xlim = c(-5,11), ylim = c(50, 62)) # plotting window

2.5 Aggregation

2.5.1 Gini Index:

Derived from a Lorenz curve. Ranges from 0 to 1, with 1 indicating that the population is uniformly distributed across surveyed rectangles, and 0 indicating that the population was recorded in one rectangle.

lordat <- lorenz_data(hlhh, yrs, qrs, species_aphia, stk_divs)
lorenz_plot(lordat) + 
  theme_minimal() +
    theme(panel.border = element_rect(colour = "black", fill = NA)) # function needs tidying

gni <- Gini(lordat)
head(gni)

# A tibble: 6 × 3
   Year Quarter `Gini Index`
  <dbl> <chr>          <dbl>
1  2000 1              0.219
2  2001 1              0.199
3  2002 1              0.231
4  2003 1              0.181
5  2004 1              0.266
6  2005 1              0.256

gni_p <- ggplot(data = gni, aes(x = Year, y = `Gini Index`)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Gini Index Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA)) 

gni_p

2.5.2 D95:

Represents the proportion of the population present in 95% of the area. Ranges from 0 to 0.95, with 0 indicating that all individuals were recorded in 5% of the surveyed area (high aggregation) and 0.95 indicating that 95% of the population were recorded in 95% of the rectangles surveyed (uniform distribution).

D <- d95(lordat)

D95 is an estimate of the proportion of the population that exists in 95% of surveyed rectangles.

d95_p <- ggplot(data = D, aes(x = Year, y = `D95`)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Proportion of Catch in 95% of ICES Rectangles in Stock Region (D95) Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA)) 

d95_p

2.5.3 Spreading Area (SA):

Derived from cumulative frequency data. Measures the area occupied by the population, taking into account variations in fish density values. High values of SA indicate homogeneous spatial distribution.

sa_data <- spreadingarea_data(hlhh, yrs, qrs, species_aphia, stk_divs)
sa <- sa_data %>%
        group_by(Year) %>%
        summarise("Spreading Area" = spreadingarea_calc(TotalNo_Dur)) %>%
        mutate(Quarter = paste(as.character(sort(unique(sa_data$Quarter))), collapse = ", ")) %>% # add quarters
        relocate(Year, Quarter) 
head(sa)

# A tibble: 6 × 3
   Year Quarter `Spreading Area`
  <int> <chr>              <dbl>
1  2000 1                   80.4
2  2001 1                   87.0
3  2002 1                   97.7
4  2003 1                   74.4
5  2004 1                   95.9
6  2005 1                   95.3

sa_p <- ggplot(data = sa, aes(x = Year, y = `Spreading Area`)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Spreading Area (SA) Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA)) 

sa_p

2.5.4 Equivalent Area (EA):

The area that would be covered by a population with homogeneously distributed density. It is equal to the mean density per individual.

ea <- sa_data %>%
        group_by(Year) %>%
        summarise("Equivalent Area" = equivalentarea(TotalNo_Dur)) %>%
        mutate(Quarter = paste(as.character(sort(unique(sa_data$Quarter))), collapse = ", ")) %>% # add quarters
        relocate(Year, Quarter)
head(ea)

# A tibble: 6 × 3
   Year Quarter `Equivalent Area`
  <int> <chr>               <dbl>
1  2000 1                    45.9
2  2001 1                    58.5
3  2002 1                    73.8
4  2003 1                    41.5
5  2004 1                    72.5
6  2005 1                    78.9

ea_p <- ggplot(data = ea, aes(x = Year, y = `Equivalent Area`)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Equivalent Area (EA) Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA)) 

ea_p

2.5.5 Spread of Participation Index (SPI):

Compares the observed spatial density distribution to the density distribution expected from a homogeneously distributed population. Ranges from 0 to 1, with 0 indicating that the population was observed in only one ICES rectangle, and 1 indicating that the population density was uniformly distributed across ICES rectangles.

SPI <- spi(hlhh, yrs, qrs, species_aphia, stk_divs)[c(1,2,4)] # function needs tidying

head(SPI)

# A tibble: 6 × 3
# Groups:   Year [6]
   Year Quarter SPI.dur
  <int> <chr>     <dbl>
1  2000 1         0.435
2  2001 1         0.450
3  2002 1         0.455
4  2003 1         0.426
5  2004 1         0.500
6  2005 1         0.483

spi_p <- ggplot(data = SPI, aes(x = Year, y = `SPI.dur`)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Spread of Participation Index (SPI) Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  theme(panel.border = element_rect(colour = "black", fill = NA)) +
  ylab("Spread of Participation Index (SPI)")

spi_p

2.6 Plot Indicator Trends

Plot all the the timeseries of indicators to visualise and compare trends.

# List of dfs
df_list <- list(cg, inert, eoo, ela, popr, poph, gni, D, sa, ea, SPI) 

# Categorise indicators
loc <- c("CoG (x)", "CoG (y)")
ran <- c("Inertia", "EOO", "ELA")
occ <- c("POPR", "POPH")
agg <- c("Gini Index", "D95", "SA", "EA", "SPI")

# Merge all outputs into a single df
sidf <- Reduce(function(x, y) merge(x, y, all=TRUE), df_list) %>%
  select(-c(nrects, nrects_p, no_haul.ids, pr_hauls)) %>% # remove some cols
  rename(EOO = convex_hull_area,
         POPR = PosAreaR,
         POPH = PosAreaH,
         ELA = `Ellipse Area`,
         SPI = SPI.dur,
         SA = `Spreading Area`,
         EA = `Equivalent Area`) %>%
  tidyr::pivot_longer(cols = 3:14, names_to = "Indicator", values_to = "Value") %>%
   mutate(Type = case_when(
     Indicator %in% loc  ~ "Location",
     Indicator %in% ran  ~ "Range",
     Indicator %in% occ  ~ "Occupancy",
     Indicator %in% agg  ~ "Aggregation")) %>%
  mutate(Indicator = factor(Indicator, levels = c(loc, ran, occ, agg)))

head(sidf)

# A tibble: 6 × 5
   Year Quarter Indicator  Value Type     
  <dbl> <chr>   <fct>      <dbl> <chr>    
1  2000 1       CoG (x)    6.35  Location 
2  2000 1       CoG (y)   55.8   Location 
3  2000 1       Inertia   14.8   Range    
4  2000 1       EOO       89.3   Range    
5  2000 1       ELA       86.2   Range    
6  2000 1       POPR       0.779 Occupancy

ggplot(data = sidf, aes(x = Year, y = Value)) +
  geom_line(aes(colour = Type)) +
  scale_x_continuous(breaks = seq(from = min(unique(sidf$Year)), to = max(unique(sidf$Year)), by = 2)) +
  facet_wrap(vars(Indicator), scales = "free") +
  theme_minimal() +
  theme(
    panel.border = element_rect(colour = "black", fill = NA),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6)
  ) +
  labs(title = "Spatial Indicator Timeseries",
       subtitle = ttl(species, stk_divs, srv, qrs, yrs)) +
  ylab("Indicator Value")

3 More Information

For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage ¹.

3.1 Author information

3.2 Acknowledgements

3.3 Software Versions

R version 4.2.2 (2022-10-31 ucrt)

FLCore: 2.6.19
ggplotFL: 2.7.0
ggplot2: 3.4.4

Compiled: Tue Jan 23 15:46:15 2024

4 References

http://flr-project.org ↩︎

Spatial Indicators Summary

Peter Kidd & Laurence Kell

23 January, 2024